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Abstract 

We present an imaging technique particularly suited to the detection of a target embedded in 
a strongly scattering medium. Classical imaging techniques based on the Born approximation fail 
in this kind of configuration because of multiply scattered echoes and aberration distortions. The 
experimental set up we consider uses an array of programmable transmitter s/receivers. A target is 
placed behind a scattering medium. The impulse responses between all array elements are measured 
and form a matrix. The core of the method is to separate the single-scattered echo of the target 
from the multiple scattering background. This is possible because of a deterministic coherence 
along the antidiagonals of the array response matrix, which is typical of single scattering. Once 
this operation is performed, target detection is achieved by applying the DORT method (French 
acronym for decomposition of the time reversal operator). Experimental results are presented in 
the case of wide-band ultrasonic waves around 3 MHz. A 125-element array is placed in front 
of a collection of randomly distributed steel rods (diameter 0.8mm). The slab thickness is three 
times the scattering mean free path. The target is a larger steel cylinder (diameter 15 mm) that 
we try to detect and localize. The quality of detection is assessed theoretically based on random 
matrix theory and is shown to be significantly better than what is obtained with classical imaging 
methods. Aside from multiple scattering, the technique is also shown to reduce the aberrations 
induced by an heterogeneous layer. 

PACS numbers: 42.25.Dd, 43.60.+d, 43.20.+g, 46.65. +g 
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I. INTRODUCTION 



Classical reflection imaging methods, such as echography or radar, are based on the same 
principle. One or several transducer (s) emit(s) a wave toward the medium to be imaged. 
The incident wave is reflected by the heterogeneities and the backscattered wave field is 
measured by the same sensor (s). The backscattered wave contains two contributions: 

• A single scattering contribution (path s in FigH]): the incident wave undergoes only 
one scattering event before coming back to the sensor (s). This is the contribution 
which is taken advantage of, because there is a direct relation between the arrival time 
t of the echo and the distance d between the sensor and the scatterer, t = 2d/c (c is 
the sound velocity). Hence an image of the medium's reflectivity can be built from 
the measured signals. 

• A multiple scattering contribution (path m in FigH]): the wave undergoes several 
scattering events before reaching the sensor. Multiple scattering is expected to take 
place when scatterers are strong and/or concentrated. In this case there is no more 
equivalence between the arrival time t and the depth of a scatterer. Thus, classical 
imaging fails when multiple scattering dominates. 

To image an heterogeneous medium, one tries to reduce the influence of multiple scat- 
tering. In that respect, multiple sensors arrays are a great improvement, since coherent 
beamforming can be achieved at emission and reception It consists in focusing the 
transmitted wave at the desired point by applying the appropriate time delays to each array 
element. In the reception mode, the same delays are applied to the received signals before 
they are summed. Single scattering signals coming from a target located at the focus add 
up coherently, whereas the summation is expected to be incoherent for multiple scattering 
signals arriving at the same time. The gain in single-to-multiple scattering provided by 
beamforming is proportional to the number of elements on the array. In medical imaging 
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where multiple scattering is usually weak at standard ultrasonic frequencies, this operation 
is generally sufficient to correctly image the medium. But in other situations, multiple 
scattering can be so high that coherent beamforming fails. The resulting echographic im- 
age is pure speckle, with no direct connection with the medium's reflectivity. There can 
be false alarms that one can wrongly attribute to the presence of a strong reflector in the 
medium. Furthermore, aberration effects distort the wave front of the focused beam, which 
may generate secondary lobes or a displacement of the focal spot. 

Our aim is to detect and image an echogene target embedded in a scatteringmedium. 



This issue has received considerable attention in the last decade 

mm 

As mentioned previously, classical imaging techniques may fail in such media because of 
multiple scattering and aberration effects. To solve this problem, various coherent interfer- 
o^etHc imaging .ech„,.. W bee„ .gg.ted B0,OlflflU. Neve.Kele.. t.ey a. 
shown to fail when the target is typically buried beneath one transport mean free path /* 
of the scattenng medium Q. Another route towards tapet detection ■„ h.gh.y scattering 
media is to tackle with the radiative transfer equation [9|, 110| . However, this approach needs 
heavy numerical computations. Moreover, the final resolution of the image is poor since it is 
limited by /* instead of half the wave length A/2. This paper proposes an original approach 
to drastically reduce the multiple scattering contribution, which can hide the echo from 
targ ets {e.g., landmines, ducts, defects) embedded in the earth [l^, in concrete structures 
(flaws, defects,...) or austenic steels [l^ for non destructive evaluation. Reduc- 
ing the influence of multiple scattering is also a challenge in optica 



(OCT) 
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18| , in seismology 19|, l20|] , in ultrasound imaging 21 



coherence tomography 



221 ] or in radar |23(]. In 



optics, correlation techniques have also been proposed to reduce the multiple scattering in 



fluence in dynamic light scattering experiments [2^, 2^, 26, 27, 28 1. However, these methods 
only address the suppression of the multiple scattering contribution in the autocorrelation 
function of intensity. On the contrary, the approach we propose here can be dedicated to a 
much wider range of applications, since it basically applies to the wave field. 

In this paper we will use ultrasonic waves in the MHz range for experimental demonstra- 
tions, but the technique can be applied to all fields of wave physics for which the multi- 
element array technology is available and provides time-resolved measurements of the am- 
plitude and the phase of the wave field. The experimental situation is the following: the 
medium we want to image is placed in front of a multi-element array (see FigH]). A pulse sig- 
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nal is sent from element i and the wave backscattered by the medium is measured by element 
j. This operation is achieved for all possible transmitter/receiver couples. The set of A^^ 
responses forms a matrix K which constitues the global response of the medium. Unlike the 
multiple scattering contribution, single scattering signals exhibit a deterministic coherence 
along the antidiagonals of the array response matrix whatever the distribution of scatterers 



291] . This particular feature can be taken advantage of to extract the single scattered waves, 



even though multiple scattering predominates. This "single scattering filter" (SSF) yields a 
filtered matrix K^, ideally devoid of multiple scattering. 

Once the separation of single- and multiple-scattered waves is performed, the detection of 



the target is achieved by the DORT method 30|, ISl] (French acronym for decomposition of 



the time reversal operator). It consists in a singular value decomposition (SVD) of the array 
response matrix. Actually the SVD is written K = UAV''", where A is a diagonal matrix 
containing the real positive singular values Aj in a decreasing order (Ai > A2 > ... > Aat). U 
and V are unitary matrices whose columns are the normalized singular vectors Ui and Vi. 
DORT has shown its efficiency in detecting and separating the responses of several scatterers 



in homogeneous or weakly heterogeneous media 



approximation and for point-like scatterers 
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311]. Indeed, under the single scattering 



331 ]. each scatterer is associated mainly with 



one significant eigenstate linked to a non zero singular value Aj. The corresponding singular 
vector Vi is an invariant of the time reversal operator KK''". Physically, each eigenvector of 
KK"!" (or singular vector of K) corresponds to a wave that, when it is sent from the array, 
focuses onto the associated scatterer. Therefore, it is possible to focus selectively on the 
corresponding scatterer and obtain its image by backpropagating Vi either physically or 
numerically. 

However, in this study, the target is hidden behind a strongly scattering slab. If we apply 
the DORT method directly to the array response matrix K, expecting that the target will be 
associated to the first singular value Ai and backpropagating numerically the corresponding 
singular vector Vi, it fails because of multiple scattering. We will show that once the 
single- and multiple-scattering contributions have been separated, DORT can be applied to 
the filtered matrix K^, and successfully detects the target despite multiple scattering. A 



detection criterion has to be applied to the first singular value Ai in order to decide i 



is detected or not. To that aim, we will refer to random matrix theory (RMT) 



to a recent work 
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a target 



351 ] and 



361 ] dealing with the statistical behavior of the matrix K in random media. 
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The efficiency of the technique will also be evaluated from RMT and shown to be better 
than classical imaging techniques. Finally, the issue of aberration will be adressed. The SSF 
is shown to strongly diminish the aberration effects which occur in scattering media. 



II. EXPERIMENTAL PROCEDURE 
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FIG. 1: Experimental setup. A 125-element array is placed in front of a random medium at a 
distance a. The whole setup is immersed in a water tank. The inter-element response kij(T,f), 
around the time of flight T and at the frequency /, is measured. It contains contributions of single 
and multiple scattering paths whose lengths belong to the interval [R — Ar/2;R + Ar/2], where 
R = cT/2 and Ar = cAt/2. Examples of a single scattering path (labelled s, grey line) and of 
a multiple-scattering path (labelled m, dashed black line) is drawn. {Xs,Zs) are the coordinates 
of the scatterer involved in path s. {Xm\ Zm) and {Xm ■, Zm) are the coordinates the first and 
last scatterers along path m. The path t represents the single scattering path associated with the 
target (continuous black line). 



The experiment takes place in a water tank. We use an N-element ultrasonic array 
(A^ = 125) with a 3 MHz central frequency and a 2.5-3.5 MHz bandwidth; each array 
element is 0.39 mm in size and the array pitch p is 0.417 mm. The sampling frequency is 20 
MHz. The array is placed in front of the medium of investigation, at a distance a = 40 mm. 
It consists of parrallel steel rods (longitudinal wave velocity = 5.9 mm//is, transverse 
wave velocity = 3 mm//is, radius 0.4 mm, p = 7.85 kg/L) randomly distributed with 
a concentration n = 12 rods/cm^. The frequency-averaged scattering mean-free path is 



7.7 ± 0.3 mm for this medium between 2.5 and 3.5 MHz 



37|. The slab thickness is L = 20 



mm. An air-filled steel cylinder with diameter 15 mm is placed behind the scattering slab. 
Our aim is to detect this echogene target. Note that the single scattered wave associated to 
the target (path t in Fig{T]) has to travel more than five scattering mean free paths through 
the random medium. Its intensity is roughly divided by exp (— 2L//e) ~ 180 as it traverses 
twice the scattering slab. Multiple scattering, in addition to aberration effects induced by 
the slab, make the detection of the target very difficult with classical imaging techniques. 
This is highlighted by the echographic image in Figj2l The first rows of scatterers in the 
slab are clearly visible. Beyond a depth of typically one mean free path (~ 5 — 10 mm), the 
image displays a speckle pattern without connection with the medium's reflectivity. The 
target, which should be visible in Figl2] around R = 70 mm, is not detected by classical 
echography. 




R [imu] 

FIG. 2: Standard echographic image of the inspected medium obtained with focused and adaptative 
beamforming both at emission and reception. The image is normalized with its maximum and 
shown in dB. 

We now turn to the acquisition of the inter-element matrix (see Fig. [1]). A 2.5- fis- 
long sinusoidal burst of frequency 3 MHz is emitted from transducer i into the scattering 
sample. The backscattered wave is recorded with the transducers of the same array. 
The operation is repeated for the A^ emitting transducers. The impulse response between 
transducers i and j is noted hij{t). An N x N response matrix H(t) whose elements are 
the A^^ impulse responses hij{t) is thus obtained. Because of reciprocity, hij{t) = hji(t) and 
H(t) is symmetric. In the following, we take as the origin of time t = 0, the instant when 
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the source emits the incident wave. 

A short-time Fourier analysis of the impulse response matrix H is achieved. The time 
signals hij{t) are truncated into successive time windows : kij{T,t) = hij(T — t)WFi(t) with 
WR{t) = 1 for t G [-At/2 , At/2], Wnit) = elsewhere. The value of At is chosen 
so that signals associated with the same scattering event (s) within the medium arrive in 



the same time window [36[. Actually, the choice of At is particularly important for single 
scattering signals, if one wants to detect scatterers properly with the DORT method. In our 
experimental configuration, we obtain a value At ^ For each value of time T, the kij 

form a matrix K. A Fourier analysis is achieved by means of a discrete Fourier transform 
(DFT) and gives a set of response matrices K(T, /) at time T and frequency /. 

III. SINGLE AND MULTIPLE SCATTERING CONTRIBUTIONS 

As an example, Figj3] shows the real part of K at the central frequency / = 3 MHz. At 
early times (Figl^a)), single scattering dominates: multiple scattered echoes have not yet 
reached the array. Figl3](b) represents K at an arrival time larger than 2(a + L)/c: at such 
times only multiple scattering can exist. K clearly exhibits a different behavior in the single 
and multiple scattering regimes. Whereas multiple scattering results in a seemingly random 
matrix K, single scattered waves exhibit a deterministic coherence along the antidiagonals 
of K. The reason for this, and its consequences on the statistical properties of the singular 
values, have been discussed in 
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36l |. We briefly recall the argument in this section. 



The signals kij{T, f) can be written as the sum of a single scattering contribution kfj{T, /) 



and a multiple scattering contribution kfUT, /) 



h,{T,f) = kf^{TJ) + kfl{T,f) (1) 

Let us express both contributions. 

The signals kfj(T, f) at a time T and frequency / correspond to the sum of partial waves 
that reach the array in the time window [T — At/2; T + At/2]. The "isochronous volume" is 
defined as the ensemble of points that contribute to the backscattered signal at a given time. 
It is formed by a superposition of ellipses whose foci are transmitter i and receiver j. In a 
far-field configuration, we can approximate the isochronous volume by a slab of thickness 
Ar = cAt, located at a distance R = cT/2 from the array and parallel to it (see Fig{T]). 
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FIG. 3: Real part of matrix K at / = 3 MHz and different arrival times T: (a) when the single 
scattering contribution is predominant (T = 58.5 /is); (b) when only multiple scattering occurs 
(T = 122.5 lis). 

For simplicity but without loss of generality, we also assume that the reflectors as well as 
the array elements are point-like. In a 2D configuration, under the paraxial approximation. 



kfAT, f) can be expressed as 



exp {j2kR) 
R 



^ Ad exp 



d=i 



jk 



2R 



exp 



jk 



[xj — Xd) 
2R 



(2) 



where k = 27r//c is the wave number in the surrounding medium, Xi = {i — N/2)p is the co- 
ordinate along the array. X^ is the transverse position of the d^^ scatterer which contributes 
to the backscattered wave at time T, the amplitude Ad depending on its reflectivity. Both 
Ad and Xd are considered random. Note that j = ^/—l in Eq|2] and has not to be mixed up 
with the subscript j which denotes the receiver index. 
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As to the multiple scattering contribution, kfj (T, /) also correspond to a sum of partial 
waves that reach the array in the time window [T — At/2; T + At/2]. They are associated 
with multiple scattering paths whose length belongs to the interval [R — Ar/2; R + Ar/2], 
where R = cT/2 and Ar = cAt/2. An example of such a path is drawn in Fig{TJ In a 2D 
configuration, under the paraxial approximation and assuming point-like transducers and 



scatterers, kfj{T, /) can be expressed as 



exp 



Jk Z, 



.(1) 



^(2) 



7(1) 7(2) 



exp 



jk 



X, 



(1) 



2Z, 



(1) 



exp 



jk- 



X, 



(2) 



2Z, 



(2) 



q=l 

(3) 

The index q denotes the g*^ path which contributes to the signal received at time T. 

are respectively the coordinates of the first and last scat- 
terers along the path q. Bg is the complex amplitude associated with path q, from the first 
scattering event at (^Xq^\ Zq^''^ until the last one at (^Xq^\ Zq^^^. 

At this stage, the theoretical expression of kfj{T,f) given in EqE] does not display any 
obvious coherence: kfj{T, f) corresponds to a sum of partial waves which are independent 
of each other since the distribution of scatterers is assumed random. One can try to express 
kfj{T, /) as a function of {xi—Xj) and (xi + xj) which corresponds to a change of coordinates 
in Eqj2j 



exp {j2kR) 
R 



exp 



jk 



[Xi 



AR 



^ Ad exp 



d=l 



jk 



2Xd) 



AR 



(4) 



deterministic term random term 

The term before the sum in EqH] does not depend on the distribution of scatterers, it is 
deterministic; on the contrary, the term on the right is random. This special feature of 
single scattering signals manifests itself as a particular coherence along the antidiagonals 
of the matrix K, as depicted by Figl3](a). Indeed, along each antidiagonal, i.e for couples 
of transmitter(z) and receiver(j) such as z -|- j is constant, the random term of EqJUis also 
constant, for any given realization of disorder. Thus, there is a deterministic phase relation 
between coefficients of K located on the same antidiagonal. It can be expressed in the 
following way : 



Pn 



ki—m.i+mi^i /) 



exp 



jk- 



R 



(5) 
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This is no longer true in the multiple scattering regime, since kfj cannot be factorized so 
simply. Note that the parabolic phase dependence along each antidiagonal of should 
be weighted by an attenuation term, decreasing with \xi — Xj\, in order to incorporate the 
directivity of transducers. Thus EqJUis not rigorously true; yet for simplicity, we will neglect 
this attenuation term in the following. 



IV. SINGLE SCATTERING FILTER (SSF) 

Now that we have explained the deterministic coherence of single scattering signals along 
the antidiagonals of the array response matrix K, we can take advantage of this special 
feature to extract the single scattering contribution from the multiple scattering background. 
Once the set of matrices K(T, /) are measured, the separation between single and multiple 
scattering contributions is achieved according to the following steps: 

• Rotation of each matrix K and construction of two sub- matrices Ai and A2. 

• Filtering of matrices Ai and A2. Two new matrices A^ and Af are obtained. 

• Construction, from A^ and A2, of the filtered matrices containing the single 
scattering signals. 

In the following subsections, we explain in details the matrix operations performed at each 
step. 



A. Matrix rotation 

A rotation of matrix elements is achieved as depicted in FigJH It consists in building two 
matrices Ai and A2 from the matrix K: 

Ai = [aiuv] of dimension (2M - 1) x (2M - 1), 

such that ai[u,v] = k[u + V — l,v — u + 2M — 1] (6) 

A2 = [a2uv] of dimension (2M - 2) x (2M - 2), 

such that a2[tt, f] = k[u + v,v — u + 2M — 1] (7) 

with M = [N + 3)/4. Here N = 125 and so M = 32 is an even number. The matrices 
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FIG. 4: Principle of the data rotation by taking the example of a matrix K of dimension N = 17. 
The black points denote the whole elements kij of K. The columns of matrices Ai and A2 contain 
the antidiagonals of K. Circles and squares represent respectively elements of Ai and A2. After 
filtering Ai and A2, the filtered matrix is obtained and corresponds to elements contained in 
the central square. Its dimension is (2M — 1) x (2M — 1) (M = 5 here). 

Ai and A2 contain the whole antidiagonals of K (see FigJlj). Therefore the coherence of 
single scattering signals now manifests itself along the columns of Ai and A2. In the next 
subsection, we will no longer make the difference between matrices Ai and A2 because they 
are filtered in the same way. They will be called indifferently A. L is the dimension of A. 
For matrix Ai, L = 2M — 1; for matrix A2, L = 2M — 2. Because of spatial reciprocity, 
K is symmetric {kij = kji). Thus, A also exhibits a symmetry: each line of its upper part 
is identical to a line of its lower part. The symmetry axis is shown as a black line in FigH] 
and corresponds to the diagonal of the matrix K. So, each column of the matrix A contains 
only M independent coefficients, even if its dimension L is superior to M. This fact will be 
crucial when the gain in signal-to-noise ratio will be assessed. 
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B. Filtering of matrix A 



The matrix A is the sum of two matrices A^ and A^, which correspond respectively to 
the single and multiple scattering contributions 



A = A^ + A 



M 



(8) 



The rotation of data can be described as the following change of coordinates {xi,Xj) — >■ 

Vu = 7^ and = ^ (9) 



V2 



V2 



In this new basis, EqJH becomes 



s exp {i2kR) 
(^uv\TJ) oc exp 



R 



■^4 



X 



(10) 



deterministic term 



^ random term 



where r„ = Zld^i^^exp 



. Each column of the matrix A^ exhibits a known 



j^' 2R 

dependence as a function of index u. On the contrary, the multiple scattering contribution 
(Eqj3]) cannot be factorized in this way. Even after rotation, the random feature remains 
along the columns and the lines of matrix A^. 

The extraction of single scattering signals can be achieved by projecting the columns of 
the matrix A on the "characteristic space" of single scattering, generated by the vector S 
whose coordinates are 



exp 



jk 



2R 



^-1/2 



The factor L ensures the normalization of S. The result P of this projection is 



P = S^A 



(11) 



(12) 



whose coordinates are 



,M 



u=l 



exp {j2kR) 



u=l u=l 
L 



M 



(13) 
(14) 



u=l 



The residual term ^^^i sJia*^ corresponds to the projection of multiple scattering signals 
on the vector S. The filtered matrix A^ is obtained by multiplying the column vector S by 
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the line vector P 

= SP = SS^A (15) 

The elements of A^ are: 

L 
m'=1 

The first term on the right-hand side of Eq^TH] is strictly equal to the single scattering 
component (Eq JTOj) . Finally, we obtain 

L 

<. = 4v + SuY.<'^u'. (17) 
u'=l 

EqJTTl can be written under a matrix formalism: 

A^= A^ + (18) 

Single scattering Residual noise 

The matrix A^ contains the single scattering contribution (A^) as wanted. But it also 

contains a residual term due to multiple scattering (SS^A^). This term persists because 

multiple scattering signals are not stricly orthogonal to the characteristic space of single 

scattering, generated by the vector S. The filtering of the single scattering contribution 

is not perfect. Nevertheless, the typical amplitude of the residual noise can be assessed. 

Since each column of A contains M independent coefficients, the filtering process decreases 

the multiple scattering contribution by a factor \fM . The single scattering contribution 

remaining unchanged, the gain in signal-to- noise ratio (in amplitude), or rather the gain in 

"single-scattering-to-multiple-scattering" ratio, is of \fM . 



exp {i2kR) 
R 



exp 



jk 



2R 



C. The filtered matrix 

Once the matrices Ai and A2 are filtered, an inter-element filtered matrix K^, of dimen- 
sion (2M — 1) X (2M — 1), is built (see FigH]) with a change of coordinates, back to the 
original system: 

• if {i — j)/2 is an integer, 

then, k^[t,j]=a([{i-j)/2 + M, {i + j)l2] 

• if {i — j)/2 is not an integer, 

then, j] = a^[{i-j - I) 12 + M, (z + j - l)/2] 
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In the following, will denote the matrix that would have been obtained if no filtering 
had been performed. K° has the same dimensions as K^, and simply contains the central 
elements of K. 

D. Illustration of the single scattering filter (SSF) 

As an example, Fig|5] illustrates the action of the SSF on experimental data. Matrices 
K''(Figj5](a)) and K^(Figl5t^b)) are shown at frequency / = 2.7 MHz and time T = 94.5 
fis. This arrival time is the one expected for the target echo. Whereas the matrix K° seems 
random, the filtered matrix displays a deterministic coherence along its antidiagonals. 
From these data, an additional operation is needed to build the image of the medium or 
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FIG. 5: Results given by the SSF at time T = 94.5 fis and frequency / = 2.7 MHz. (a) Real part 
of matrix K''. (b) Real part of filtered matrix K^. 
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rather, to detect and image the target placed behind the scattering slab. To that aim, several 
imaging techniques are compared in the next section. 

V. DETECTION AND IMAGING OF THE TARGET 

In this section, we apply the filtering procedure described above to the detection and 
imaging of a target embedded in a scattering medium. At a given frequency, two imaging 
techniques are compared: focused beamforming (equivalent to echography in the frequency 
domain) and the DORT method. As we will see, there is no interest in combining the 
SSF with focused beamforming (FB). But its combination with the DORT method provides 
excellent results. 

A. Focused beamforming (FB) 

The simplest way to image the target is to achieve a direct backpropagation of the mea- 
sured signals K(T, /), for a given time and frequency couple. The focal plane is parallel to 
the array and located at depth R = cT/2; it is discretized in a set of points. The back- 
propagation algorithm is based on the Born approximation. It consists first in calculating 
the propagation operator G, whose elements are the spatial Green functions gn between the 
^th g^j-j-Qy element and the point in the focal plane, as shown in FiglHl The medium is 
considered as homogeneous with a wave celerity c equal to that of the surrounding medium. 
At a given time T (corresponding to depth R = cT/2) and frequency /, the final image is 




R=cT/2 



FIG. 6: Principle of FB. The focal plane is parallel to the array and located at depth R = cT/2 
from it. It is discretized with a sampling period equal to the array pitch. The distances rn are 
much larger than the wavelength. 
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a vector I, the absolute value of the backpropagated wave field, which can be plotted as a 
function of X, the transverse coordinate in the focal plane: 

I=|GtK°G*| (19) 

This backpropagation algorithm is the equivalent of echography in the frequency domain, 
with a poorer temporal resolution due to the duration At of the time-windows. In the 
following, we will refer to this imaging technique as "focused beamforming" (FB). FB has 
the same drawbacks as classical echography (Fig|2]), particularly the presence of speckle 
which hides the target. We can point out its inability to detect the target by considering 
the image obtained at time T = 94.5 /is and / = 2.7 MHz (see Figl?]). This arrival time 
corresponds to the target depth, and 2.7 MHz is the frequency for which the mean-free path 
of .he Slab the .a.,est .uMp.e seaUe.,, is the weakest, Q. The presence of the 
scattering slab seriously degrades the image. Two peaks seem to arise but neither of them 
is located at the expected position. 




A" [imn] 

FIG. 7: Image obtained by FB at time T = 94.5 and frequency / = 2.7 MHz (black line). The 
ideal image obtained without the scattering slab is in grey. The vertical black line indicates the 
position of the target. 

The resulting image I (Eq JT9l) could also be averaged over the whole frequency domain 
(the result would be comparable to Figl2]), or on specific frequency bands for which the 
detection is more likely. Indeed, it is possible to establish a detection criterion based on 
speckle statistics, for a given probability of false alarm. This will be done, as well as for 
other techniques, in SecJVIl 
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In order to improve the results provided by FB, one could think of applying Eq{T9] to 
the filtered matrix instead of the raw matrix K°. Yet it can be shown (Appendix 
that this would not change the result. A short interpretation can be given. FB relies on the 
fact that single scattering signals will add up coherently as long as they come from a focal 
point at depth R = cT/2 (axial focusing) and the desired transverse position X (lateral 
focusing). The SSF also enhances single scattering signals associated with scatterers located 
around R = cT/2, but independently from their transverse position X. Now, if we build the 
echographic image (Eq lTQl) from the filtered signals, there is a redundancy in the choice of 
depth R; consequently, the SSF does not bring anything when it is followed by FB. The two 
techniques are not complementary. A more rigorous demonstration is given in Appendix \M 



B. DORT applied to 



DORT 30, ISjJ consists in achieving the singular value decomposition (SVD) of the array 



response matrix before imaging the medium: 

K° = U°A°V°t (20) 

where A° is a diagonal matrix containing the real positive singular values A° in a decreasing 
order (A? > A2 > ... > A^j). U° and V° are unitary matrices containing the normalized 
singular vectors U? and V?. Under the single scattering approximation, each scatterer of 
the medium is mainly linked with one singular space associated to a non-zero singular value 
A°. The corresponding singular vector V? is the signal to apply to the array in order to 
focus on the corresponding scatterer. Thus, the numerical backpropagation of the singular 
vectors allows to image each detected scatterer. The image provided by DORT is a vector 
I?(T, /) which corresponds to the absolute value of the backpropagated wave field 

I0 = A°|V?G*| (21) 

I?(T, /) represents the backpropagated image at time T and frequency / of the z*'^ singular 
vector VP(T,/). 

We would like the first singular space (linked to A^) to be associated with the target 
echo. However, the forest of rods in front of the target results in multiple scattering which 
hides the target echo. Its influence is illustrated in FigJHl which displays the result obtained 
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with the DORT method at the expected arrival time for the target echo (94.5 /is), and at 
the frequency for which multiple scattering is at its weakest (2.7 MHz). Figj8](a) represents 

(a) 
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FIG. 8: (a) Real part of the first singular space A^U^V^^ of K° obtained at time T = 94.5 ^s and 
frequency / = 2.7 MHz. (b) Image obtained by backpropagation of the first singular vector at 
the same time-frequency couple. The DORT image (black line) is compared with the ideal image 
obtained without the forest of rods (grey line). The vertical black line indicates the position of the 
target. 

the real part of the first singular space A^U^V"^ of the matrix K° (see FigJS](a)). It does 
not display the feature of a single scattered echo [i.e concentric circles centered around 
the target position like in FiglHl^a)). Multiple scattering results in a random matrix K'' 
(see Figj5](a)) whose singular spaces are random, without connection with the direct echoes 
of scatterers located in the isochronous volume. The corresponding image obtained by 

18 



backpropagation of the singular vector V° is shown in Fig|8]^b). No peak is observed around 
the target location. 

C. DORT applied to the filtered matrix 

Here we combine the DORT method with the SSF described in Sec JIVI (SSF+DORT 
approach). The procedure is the same as the one described in the previous subsection, 
except that is replaced by K^. Figj9] illustrates the success of this combination. The 




A' [iimi] 

FIG. 9: (a) Real part of the first singular space U^V^^ of obtained at time T = 94, 5 fis 
and frequency / = 2, 7 MHz. (b) Image obtained by backpropagation of the first singular vector 
Vf at the same time-frequency couple. The image (black line) is compared with the ideal image 
obtained without the forest of rods (grey line). The vertical black line indicates the position of the 
target. 
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filtered matrix (see Figl5]^b)) already shows a possible feature of a single scattered echo 
coming from the target. Nevertheless, it is still perturbed by a residual multiple scattering 
contribution. Once the SVD of is achieved (K^ = U^A^V^''^), the first singular space 
Af U^Vf^, shown in FiglHI^a), clearly exhibits the feature of the single scattered echo coming 
from the target. The backpropagation of the singular vector is shown in Figin](b). 
The image clearly displays a peak at the target position, with a spatial resolution that is 
comparable to the free-space situation! Yet the peak amplitude is lower, since the intensity 
of the coherent wave coming from the target has undergone an attenuation of ~ exp(2L//e) 
due to the forest of rods. 

The comparison of Figsd [S] & [3 illustrates the success of our approach. The SSF elim- 
inates a major part of the multiple scattering contribution. It allows the DORT method 
to extract properly the target echo, which was not possible with classical imaging tech- 
niques. However, these results are obtained for a time-frequency couple chosen arbitrarily. 
An imaging procedure must work blindly, without knowing the depth of the target or the 
frequency band in which the medium is less scattering. The detection of the target has to 
be systematized in order to select automatically the arrival time(s) and frequency band(s) 
for which the target is detected and can be imaged. Moreover, multiple scattering signals 
can generate false alarms because of speckle fluctuations that one can wrongly attribute to 
the presence of a strong reflector in the medium. Thus, whatever the imaging procedure, 
a rigorous detection criterion has to be established in order to discriminate artifacts and 
compare the different techniques on a common basis. This is done in the next section. 

VI. DETECTION CRITERIA 

At a given time T and frequency /, a target will be detected if the observed quantity is 
above a certain threshold. In the case of the DORT method, the detection criterion will be 
applied to the first singular value Ai. In the case of FB, the relevant variable is the maximum 
of the image I (Eq lTOj) . Since the scattering medium is considered as one realisation of a 
random process, setting the detection criteria requires a statistical model for the probability 
density function of Ai and of I. Then a probability of false alarm (PFA) is fixed, and the 
corresponding detection thresholds can be established for both methods, which allows to 
compare their results for the same PFA. 
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The statistical behavior of the singular values Aj and of the echographic image I in the 
multiple scattering regime has to be known. To that aim, we have performed the same 
kind of experiments as described in FiglH except that the target has been removed. The 
experimental procedure remains unchanged and a set of matrices K°(T, /) and K^(T, /) is 
obtained. 

We first consider the DORT method, and the statistical properties of the singular values 



iC^ in connection with random matrix theory (RMT), as discussed in recent 



of and 
papers Q 



36| . Experimentally, before achieving the SVD, the matrices K° and of 
size (2M — 1) X (2M — 1) are truncated into matrices of size M x M by keeping only one 
element in two. This operation is needed in order to remove short-range correlations that 



may exist between adjacent entries 36|]. These correlations have an important influence on 
the distribution of singular values and removing them simplifies the problem. For the sake 
of simplicity, we will continue to note the truncated matrices and K^. The SVD of these 
matrices is achieved and a set of M singular values A°(T, /) and Af (T, /) is obtained at each 
time-frequency couple (T, /) . The singular values are normalized by their quadratic mean : 

A. = , (22) 

/ 1 S-^I \2 

This normalization allows to meet the hypothesis usually made in RMT which consists in 



assuming a variance of for matrix coefficients 36| . Once this normalization is performed 
at each time-frequency couple, a histogram of dimensionless singular values is obtained by 
averaging over time T, frequency / and rank i. Two estimators, po(A) and Pf(A), of the 
singular values distribution are finally obtained, respectively for matrices K° and K^. The 
results are displayed in FigJTUl 

In the multiple scattering regime, once short-range correlations are removed, we expect 
the matrix K° to be random. In that case, RMT predicts t hat the distribution of singular 
values should follow the "quarter-circle law" (for M >> 1) 
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, ^VI^ for < A < 2 
PQc(A) = <; - (23) 
elsewhere 

As pointed out by FiglTOlfa). the experimental distribution of singular values Pn(X) deviates 



from the quarter circle law pQc{X). The reasons for that have been discussed in 3^. When 
a detection threshold is fixed for the first singular value A?, we will use experimental data 
Po rather than the theoretical quarter-circle law (pqc)- 
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FIG. 10: (a) po{X) (white disks) is compared with the quarter-circle law pQc(A) (black line, Eq l23p . 
(b) /5_f(A) (white disks) is compared with the Hankel law phW (black line), (c) pi{X) (white disks) 
is compared with the Rayleigh law pr{X) (black line, Eql25p. 
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Contrary to K°, the filtered matrix is characterized by a deterministic phase relation 



along its antidiagonals. This kind of matrix has already been studied [36|]. displays the 
same statistical properties as a Hankel random matrix. A Hankel matrix is a square matrix 
whose elements belon ging to the same antidiagonal {i + j = constant) are equal. In the 



literature, Bryc et al. (39| have proved, for normalized random Hankel matrices, the almost 
sure weak convergence of the distribution of singular values to a universal distribution of 
unbounded support p_ff(A). In the following, the distribution will be referred to as 

the "Hankel law". To our knowledge, no analytical expression of the Hankel law has ever 
been found and only a numerical simulation can provide an estimate of Ph{^)- In FigJTOlfb). 
the experimental distribution of singular values of K^, p_f(A), is compared to the Hankel 
law. The agreement between both curves is excellent. Thus, we will rely on the statistical 
behavior of Hankel random matrix, when a detection criterion is set on the first singular 
value Af . 

We now consider the echographic image I(T, /) and build an estimator for its probability 
density function. Experimentally, the points where the reflectivity of the medium is esti- 
mated have to be chosen carefully. The image vector I(T, /) has to display independent 
coordinates. So, each one has to be associated with a different resolution cell. In pratice, 
we have considered the points located at the same transverse position as the transducers 
implied in the truncated matrices K° and K^. The points of the image are separated by 
a distance 2p ~ 0.84 mm larger than the size of the resolution cell ^ ~ 0.76 mm (with 
D = (2M — l)p the array size). Each image I(T, /) contains M independent coordinates 
Ii{T,f). These coordinates are normalized by their quadratic mean at each time T and 
frequency /: 

h{TJ) = ^J£l£= (24) 

Once this renormalization is performed, a histogram of the dimensionless image can be built, 
averaging over all time- frequency couples. The estimator p/(A) of the image probability 
density function is plotted in Fig JTUT c). 

In the multiple scattering regime, we expect // to be the modulus of a gaussian complex 



random variable with zero mean and variance unity 40|, |41|. The associated density of 
probability is the Rayleigh law p_r(A): 

PR{X) = 2Xexp{-X^) (25) 
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Pr{X) is compared to the experimental estimator p/(A) in FigJTOl The agreement between 
theory and experiment is excellent. Consequently, the Rayleigh law Pr{X) will be considered 
when a detection criterion is set for the main peak of the echographic image. 

Now that the probability density functions of A^, Af and /; are known, a detection crite- 
rion can be set for each imaging technique. The relevant quantity for that is the distribution 
functions Fi of the first singular values, Af and A?, and of the main peak of the echographic 
image. Fi will directly provide the probability of false alarm PFA for the target detection 
issue, since PFA{a) = 1 — Fi{a) = Prob {a < A} (where a is the detection threshold and 
A is the variable on which the detection threshold is applied). 

As we have seen, the distribution po{X) does not strictly follow the quarter-circle law. 
Thus, the distribution function F^ of A^ will be estimated from experimental measurements 
rather than from an analytical expression. By building the histogram of the first singular 
value A^ and then considering its primitive, we obtain an estimator of the distribution 
function F^ of Xi, which is plotted in FigJTTl 

In the case of K^, the distribution of its singular values was found in good agreement 
with the Hankel law. The distribution F^ of the first singular value of a Hankel random 
matrix is calculated numerically jsol. F^{X) is also plotted in FiglTTl 

As to FB, the Rayleigh law prQCi was found to fit the data properly. Unlike the sin- 
gular values of a random matrix 
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431], the coordinates of the echographic image are 
independently distributed. In that case, the distribution function Ff-{X) of the main peak. 



max = max 



of the image can be directly deduced from Pr{X). F{^{X) is equal to the 
M^^ power of the distribution function F^{X) of any image coordinate: 

FfiX) = Prob {/max < a} = [F^iX)f , 

with, F^(A) = [ dxpR{x). 
Jo 

The distribution function -F/^(A) is plotted in FigJTTl 

At this stage, we have reliable models for the distribution function Fi of the relevant 
variable for all three techniques (DORT method, DORT method combined with the SSF, 
FB). An admitted probability of error 7 (z.e, a false alarm rate) is chosen. The three 



detection thresholds a can be obtained from 



36j: 



a = Ffi(l-7) (26) 
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FIG. 11: Distribution functions -F°(A) (dash-dotted line), F(^{X) (dashed line) and Ff {X) (contin- 
uous line) obtained for M = 32. The vertical lines represent the detection thresholds a for a PFA 
7 = 10-3. 

TABLE I: Table of detection thresholds deduced from Eql2Sl taking 7 = 10^^ and of detection 
conditions established in Appendix [Bl with M = 32. 



Imaging technique DORT, K° DORT, FB 

Detection threshold a = 2.39 a = 2.69 a = 3.15 

Detection condition ^ > ^ 0.42 ^ > ~ 0.12 > ^ ~ 0.14 



In FigJTTl the detection thresholds are represented with vertical lines; the admitted PFA 7 
has been set to 10"^ for all three imaging techniques. The corresponding numerical values 
are given in Tab HI 

Once the detection thresholds are known, we can also evaluate the performances of each 
technique for detecting a target. It consists in predicting the signal-to- noise ratio above 
which the target is detected ("noise" meaning here multiple scattering). Let cr|. and cr|,j be 
the power of signals associated with the target and the multiple scattering contribution. We 
can predict above which ratio the target is detected by each technique with the same 
probability of false alarm 7 = 10^'^. Details of calculations are given in Appendix [Bl The 
performances of each technique are summarized in Tab HI 
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The SSF+DORT approach is the most efficient in terms of detection. Its detection hmit 
decreases in as for FB, which is better by far than the classical DORT method whose 
detection limit decreases in For a given probability of false alarm, the SSF+DORT 

approach succeeds in detecting the target for slightly smaller signal-to-noise ratios than FB 
(0.12 vs 0.14). As we will see in Sec lVIIl this is not the only reason why this approach 
provides better results than FB: it also diminishes aberration distortions, whereas the image 
obtained with standard echography can be strongly damaged by aberration. 



faj (b) 




f[MHz] f[MHz] 

FIG. 12: Time-frequency evolution oi Imax{^o\)) , A]* (middle), Af^(bottom). On the left, no criterion 
of detection is applied. On the right, the gray scales have been adjusted so that for each technique, 
all values below the detection criterion are represented in black. The probability of false alarm is 
the same for all techniques (7 = 0.1%). 

The detection thresholds summarized in TabUare now applied to the experimental results. 
Figll2] represents the time-frequency evolution of Imax-, A? and Af . From the figures on the 
left, it is difficult to decide for which time T and frequency /, the target is detected. The 
application of the detection thresholds, based on the same probability of false alarm, provides 
an unambiguous answer (see Fig JT2|) . The target is detected over very few time- frequency 
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couples for FB and the classical DORT method. On the contrary, the combination of DORT 
with the SSF manages to detect the target over the frequency band 2.65 — 2.8 MHz and a 
7-/is-long temporal window. It is no accident that the target is best detected around 2.7 
MHz. Actually, the forest of rods exhibits a larger scattering mean free path /g (~ 10 mm) 
around this frequency 37|. Thus, the slab is more transparent in this frequency bandwidth, 
the direct echo of the target is less attenuated by scattering. Setting detection thresholds 
with the same PFA provides a systematic way to compare the three techniques and detect 
the target, which would not have been possible by a simple look at the echographic image 
(Figl2]) or at singular values (FigHSl left column). 
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FIG. 13: Images of the target: (a) reference image, without the forest of rods; (b) FB through 
the forest of rods; (c) DORT combined with the SSF through the forest of rods, (d) DORT alone 
through the forest of rods. All images have been renormalized by their maximum. 



Now that we have determined the time-frequency couples for which the target is detected, 
the final image can be obtained by summing, for each time T, the images over the frequencies 
/ which fulfill the detection criterion. The final image is displayed as a function of the 
transverse position X and the depth R = cT/2. FiglTSl shows the images obtained for each 
imaging technique. The "ideal" image obtained without the forest of rods is also shown 
and constitutes the reference (Fig JTST a)). The results are excellent (see Fig JTST c)): the SSF 
provides an image of the target which is comparable to the reference image, although its 
axial resolution is a bit degraded: the temporal spreading of the target echo compared to the 
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"ideal" case is due to the loss of a major part of the initial frequency bandwidth. Whereas 
the emitted signal displays a frequency bandwidth of 1 MHz, the target echo is only detected 
over a bandwidth of 0.15 MHz. The lateral resolution is nearly as good as the "ideal" image 
(see FigJTWb)). and the correct position of the target is obtained. 

The other techniques (DORT alone, FB) manage to detect the target, but only in very 
narrow frequency bands, and they are strongly affected by aberration. Even though one 
peak is observed around the expected location of the target, there are secondary lobes (for 
both techniques) and a displacement of the focal spot (for FB) (see Fig jMT a)). This is due 
to the inadequacy of the Born approximation: when backpropagating data, the medium 
is considered as homogeneous, which is obviously not valid here. The various techniques 



of aberration corrections 
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471] are difficult to apply in our experimental con- 



figuration, because of multiple scattering. On the contrary, the SSF+DORT approach is 
less sensitive to the error made when backpropagating data in a supposedly homogeneous 
medium. As we will see in the next section, the filtering of the antidiagonals of K smooths 
the distortions endured by the wave front, which tends to diminish the aberration effects. 

An other interesting observation is the occurrence of an echo above the detection threshold 
around T = 115 /is (see FigJTWf)). This echo is not an artifact due to multiple scattering. 
It corresponds to circumferential waves that have propagated around the air-filled cylinder. 
This phenomenon has been already observed with the DORT method 48|]. The difference of 
arrival times between the specular echo (90 fis < T <97 /is) and this second echo (T ~ 115 
/is) is compatible with an Aq Lamb mode. This is interesting for a better characterization 
of the target. 



VII. ABERRATION 



So far, we have dealt with the issue of multiple scattering, and how it could be partially 
eliminated by a matrix manipulation, in order to improve target detection. But multiple 
scattering is not the only enemy in imaging and detection. Even if only single scattering 
takes place, a heterogeneous layer such as the forest of rods induces aberrations that distort 
the wavefront reflected by the target. To reconstruct an image, both FB and DORT rely on 
the hypothesis that the medium has a constant speed of sound, which is clearly not true. 
In this section, we examine the impact of the SSF on aberrations. We consider the same 
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FIG. 14: Section of the images obtained at the depth of the target {z = 68 mm). The reference 
image without the forest of rods is in continuous black Hne. (a) DORT method alone (continuous 
grey line) and FB (dashed line), (b) DORT method combined with the SSF (continuous grey line). 
All curves have been renormalized with their maximum. 

experimental set-up, except that now the aberration undergone by the target echo will be 
examined independently from the multiple scattering contribution of the forest of rods. To 
that end, the impulse response matrix has been measured in three configurations: 

• Configuration 1: With the scattering slab and the target; it corresponds to the exper- 
imental situation studied until now. 

• Configuration 2: With the scattering slab alone (the target has been removed). 

• Configuration 3: With the target alone (the scattering slab has been removed). 
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Let H(')(t), denote the corresponding matrices, where the superscript i stands for the 
configuration (1, 2 or 3). In order to investigate aberration eff'ects apart form multiple 
scattering, we calculate the matrix H = H*^^-' — H*^^-'. H contains only signals linked to 
the target. Particularly, its first arrivals correspond to the single scattering contribution 
(ballistic) coming from the target. The later echoes correspond to multiple scattering paths 
involving both the target and the forest of rods. In FigJTSl a line of the matrix H(t) is 
displayed as a function of the arrival time t in the time window 90 — 100 /is. It is compared 
with the same line of matrix tl^^\t). The comparison of matrices H and H*^^) highlights the 
phase and amplitude distortions undergone by the wave-front refiected by the target when 
it comes through the scattering slab. 




FIG. 15: (a) Line 64 of matrix H(^)(t) in the time-window 90 — 100 /zs. (b) Line 64 of matrix il{t) 
in the same time-window. 
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Let us study the action of the SSF on the distortions. The frequency spectrum of H and 
H^^^ is calculated by means of a DFT in the time-window 90 — 100 fis. Two matrices, K(/) 
and K'^^)(/), are obtained at each frequency /. The distortions induced by the scattering 
slab can be quantified by a matrix D(/) whose coefficients dij{f) are 

If the scattering slab had no effect, the coeffients dij would be real (no phase distortion) and 
equal to unity (no amplitude aberration). Obviously, in our case, the coefficients dij show 
both amplitude and phase distortions (FiglTGll. 




J 

FIG. 16: (a) Real part of the distortion matrix D'' at the frequency / = 3.1 MHz. (b) Real part 
of the filtered distortion matrix at the same frequency. 
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Let us express the coefficients kij and klj\ Considering Eqj2]with only one scatterer (the 



target itself) we have: 



kifif) exp 



jk 



2Rr 



exp 



2Rr 



(28) 



where the coordinates {X^, Rt) correspond to the target location. For the sake of simplicity, 
we have removed the phase term '^^p(i^^^T) ^^^q reflectivity term in Eq|2], which are 
unimportant here. Using Eql27] and Eq j28l we have 



hjif) oc dij{f)exp 



jk- 
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2R 



(29) 



As seen in Sec lIVl the SSF consists in projecting the antidiagonals of the matrix K on the 
characteristic space of the single scattering contribution. Once this operation is performed, 
the coefficients of the filtered matrix can be expressed as (see Appendix [C]) 



oc ei+^„i(/)exp 



where the coefficients are given by 



if V is an odd number, 
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then, e„ 
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(32) 



and the symbol < . > denotes an average over the variable in the subscript. Therefore 
the coefficients e„ result from a smoothing of dij: the effect of the SSF is to average the 
distortion coefficients along each antidiagonal. Fig JTGT b) represents the ffitered distortion 
matrix, D^. Its coefficients are dfj = ej+j_i. The comparison of matrices D° and shows 
that the SSF reduces the fluctuations of the distortion coefficients (see Fig JTGj) . 

The standard deviation std [cp] of coefficients Cp is smaller than std [dij] by a factor of 



\/Nind- Nind is the number of independent elements along each antidiagonal of matrix D, 
which depends on the coherence length of the aberrator relatively to the array pitch. If the 
matrix D contained independent entries, Nind would be equal to M. In our experimental 
configuration, there are correlations between neighbouring elements belonging to the same 
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antidiagonal, and Nind ~ "y- Fig{T7] illustrates the action of the SSF on the fluctuations of 
the distortion coefficients. The evolution of the ratios \{dij) \ /std [dij] and |(ep)| /std [cp] is 
plotted as a function of frequency. Note that |(ep)| = \ {dij)\. The ratio between both curves 
is close to y/Nind — \J^- = 4. 

1.5 
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FIG. 17: Ratios /std [djj] (continous grey line) and |(ep)| /std [cp] (black continous line) are 

plotted as a function of frequency. The horizontal dashed line corresponds to the detection threshold 
of EqIMl 

Now that it is clear that the SSF has a smoothing effect on the aberrations, let us evaluate 
by how much it will improve target detection. We start by keeping only one in four elements 
of the filtered matrix K^: its new dimensions are ^ x As before (see Sec lVip . this is done 
to remove short-range correlations between matrix elements, which allows us to use relatively 
simple results derived from RMT. Next, we perform a singular value decomposition: 

= U^A^V^t (33) 

If there were no aberrations at all, backpropagating as usual the first singular vector 
would focus at the target position. This is not the case here, because of aberrations induced 
by the scattering layer. In order to analyse their effect, let us write as the sum of a 
"smoothed" matrix (K^) and a perturbation AK^: 

K^= (K^) + AK^ (34) 

Matri^ToFrank 1 Random Hankel matrix 
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Note that the absence of subscript behind the symbol < . > means that we now consider 
ensemble averages. From Eq j30l we have the coefficients of (K^) and AK^: 
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2Rt 



exp 
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(35) 



(36) 



Note that (K^) may be written as 



(37) 



with (cp) the average distortion and K^^^ the response of the target without the scattering 
slab. (K^) is of rank 1 and its only singular vector focuses at the exact location of the 
target. AK^ is the perturbation due to the aberration effects: it explicitly depends on the 
fluctuations of the distortion coefficients Cp around their mean (cp). AK^ has a particular 
feature: since the coefficients [ei+m-i(/) — {^i+m-iW))] constant along each antidiagonal 
{I + m =constant), there is a deterministic phase relation between the coefficients of AK^ 
located on the same antidiagonal. As shown in 36|], this kind of matrix has the same 
statistical behavior as a random Hankel matrix. Backpropagating the ffist singular vector 
of will focus at the target position, as long as the perturbation AK^ is weak compared 
to (K^). a threshold can be evaluated, based on RMT (the details are given in Appendix 
[Pll . The SVD will successfully extract (K^) as long as 



Kep)| 



> 



a 



(38) 



std [cp] ^M/2 

with a the detection threshold found for the ffist normalized singular value in the case of a 
random Hankel matrix (see Sec I VII) . This detection threshold is displayed with a horizontal 
line in FigJTTl a has been calculated here considering the distribution function F^{X) 



obtained for a random Hankel matrix of size ^ ^ "^i^h ^ =16. As before, the PFA 
has been fixed to 7 = 10"'^ and we have obtained numerically the threshold a = 2.52. The 
SSF+DORT approach succeeds at the frequencies / where the ratio | (cp) | /std [cp] (black 
curve) is above the detection threshold defined in Eqj38] (see Fig JT7|) . 

As an illustration, we choose to work at / = 3.1 MHz, a frequency for which we are 
just above the threshold of EqlHHl The SVD is applied to the matrices K° and K^. The 
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unwrapped phases of the first singular vectors V° and are plotted in Fig JTST a). They 
are compared with the "ideal" phase which is obtained without the scattering slab. It is 
given by the parabolic term k ^^'2Rt^ which allows to focus on the target when the first 
singular vector is backpropagated numerically. The unwrapped phase of the first singular 
vector is a relevant observable because it directly shows the phase distortions of the wave 
front which focuses on the target. Without prior filtering, the strong phase distortions of 
the wave-front result in a first singular vector V° whose unwrapped phase exhibits erratic 
deviations from the "ideal" parabolic law. On the contrary, the SSF leads to a first singular 
vector whose unwrapped phase is very close to the "ideal" case (FiglTSlfa)). 

Fig JTST b) displays the images obtained with the numerical backpropagation of the singular 
vectors on the focal plane. The backpropagation of Vf focuses nicely on the target, nearly as 
if there was no aberrating slab. On the contrary, without prior filtering, the singular vector 
V° does not focus on the target and it is impossible to deduce the target location from 
the image. Finally, we also show the result of FB applied to the matrix K° (see Sec JV A[) . 
Once again, the phase distortions induced by aberration are so large that the image displays 
several main lobes with no connection with the target location. 

In this example, the action of the SSF on phase distortions is obvious. The fluctuations 
of the distortion coefficients dij are diminished by a factor of -^Z^- As long the criterion 
(Eq j38l) is fulfilled, the SVD succeeds in extracting the unaberrated part (K^) from the 
measured matrix. Then, backpropagting the first singular vector provides the correct 
target location as shown in FigJTSlfb). 

In the previous sections, we have compared the SSF to FB and showed that it gave much 
better results when trying to detect a target behind a multiple scattering and aberrating 
layer. The smoothing of the wave-front distortions provided by the SSF partly accounts for 
this. Even if the single scattered echo coming from the target is sufficiently large compared 
to the multiple scattering contribution, FB may fail in detecting the target because of 
the strong phase distortions induced by the scattering slab. On the contrary, the SSF 
smoothes the phase distortions enough for the SVD to extract the undistorted wave-front. 
Backpropagating the first singular vector in a virtually homogeneous medium finally allows 
to image the target, with no secondary lobes, and no displacement of the focal spot. 
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FIG. 18: Influence of the phase distortions on the three imaging techniques at frequency / = 3.1 
MHz (a) Unwrapped phase of the first singular vectors (black dashed line) and (gray dashed 
line) compared with the "ideal" wave-front (continuous black line) which would be obtained if the 
slab was removed, (b) Images of the target obtained with FB (black dashed line), DORT applied to 
(grey dashed line) and DORT applied to (grey continuous line). The ideal image obtained 
without the scattering slab is in continuous black line. All curves have been renormalized with 
their mean amplitude. 

VIII. CONCLUSION 

The approach we developed here combines a "single scattering filter" with the DORT 
method. It greatly improves the performance of an array of transmitters/receivers in de- 
tecting and imaging a target hidden behind a scattering medium. On the one hand, the 
short time-frequency analysis allows to select the frequency bandwidth (s) favourable to the 
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detection of the target, unlike classical echography which is performed in the temporal do- 
main. This has been made possible by setting a detection criterion based on random matrix 
theory. We showed that removing most of the multiple scattering contribution significantly 
improves the performance of the DORT method in random scattering media. The results 
are even better than focused beamforming in terms of target detection. On the other hand, 
the SSF+DORT approach is shown to strongly diminish the influence of aberration effects 
(secondary lobes, displacement of the focal spot) which burden classical imaging techniques. 
The perspectives of this study are numerous. This technique can be applied to other types of 
waves (electromagnetic, seismic, etc.) as long as a coherent array of independent elements is 
available. A future step will be to test this approach in real situations, such as the detection 
of a target embedded in the soil or of cracks in concrete structures, steel blocks etc. 



IX. ACKNOWLEDGMENTS 



The authors wish to acknowledge the groupe de recherches IMCODE of CNRS (GDR 
2253). 



APPENDIX A: FB APPLIED TO 

The aim of this appendix is to show why FB and the SSF are not complementary. 
Using EqUni and under the paraxial approximation, the coordinate /; of the echographic 
image at the arrival time T = 2R/c and frequency / can be expressed as 



i 3 ^ ^ 



(Al) 



One can write Eq JAll in the basis defined in EqlH 

k 



Y Yl ""'^ ^^p { ~^jr + (^^ ~ } 



(A2) 



where the coefficients a^v are deduced from kij after the data rotation described in Sec JIV A] 
Eq iA2l can be simplified into 



/, {T,f) = \VLY exp { -J ^ [ {y. - v^^O '] } E < 



(A3) 



Pv 
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where the coordinates Su of the vector S are given by EqJTTl We see here that the principle 
of FB can be decomposed into two steps corresponding respectively to the sum over index 
V and the sum over index u. Actually, the second sum corresponds to the projection of the 
columns of matrix A on the characteristic space of single scattering generated by the vector 
S, as seen in Sec JlV B] Thus, the SSF constitutes one of the two steps of FB. 

Now, one can try to combine the SSF with FB. We can express the image lf{T, f) that 
we would obtain from the filtered matrix K^. It can be deduced from Eg J A3 1 replacing 
by the filtered elements a^^, 



lf{T,f) 



k 



V2: 



(A4) 



Using the expression of given in Eq{T7l If {T, f) becomes: 



ifiTJ) 



H ^ <- -I ' T» 



+ J^exp <^ -J 



2R 



V-V2: 



Xi 



u'=l 



As the vector S is normalized, we have s* s^j = 1 and the latter equation can be simplified 
into 



expi-J^ 



k 



(^yy - a/2: 
- \p2xi 



u 

E^«"^ 



This equation is strictly identical to Eq lA3[ 

Ii{TJ) = lf{TJ) 



(A5) 
(A6) 

(A7) 



which means that the images built from the raw data (matrix K) or the filtered data (matrix 
K^) are identical. There is no interest in combining the SSF with FB. 



APPENDIX B: DETECTION CONDITION IN PRESENCE OF MULTIPLE 
SCATTERING 

The aim of this appendix is to predict the performances of each technique (FB, DORT 
alone, DORT combined with the SSF) in detecting a target hidden behind a diffusive slab. 
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Let erf. and afj be the power of signals associated with the target and the multiple scattering 
contribution. The performance of each technique is assessed by determining the signal-to- 
noise ratio — above which the target will be detected. 

We consider a time of flight T corresponding to the arrival time for the target echo. The 
measured matrix K° (of dimension M x M) can be decomposed as: 



K° = + K 



M 



(Bl) 



j^M corresponds to the multiple scattering contribution; its coefficients are assumed to be 
gaussian complex random variables, identically and independently distributed, of variance 
a1,j and with zero mean. Because of spatial reciprocity, K'^ is symmetric. K'^ is associated 
with the target echo. Its coefficients can be expressed as : 



aT exp 



jk 



[Xi 



2Rn 



exp 



jk 



{xj - Xt) 
2Rt 



(B2) 



where {Xt, Rt) are the coordinates of the target. The aberration effects generated by the 
diffusive slab are neglected. 

1. Focused beamforming 



The coordinates of the image I are given by 



(B3) 
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Let us express the intensity of the image, using the decomposition /c° = kfj + kf-: 
If = kpgkj* exp |-^^ [i^P - ^if - (^r - xif + {xs - xif - {xt - xif] | (B4) 



p,q,r,s 



^pg^^s* exp |-^^ [i^p - ^if - i^r - xif + {xs - xif - {xt - Xif] | (B5) 



p,q,r,s 



+ ^n^rl* exp [(^p - ^if - (^r- - xif + (x, - xif - {xt - xif] | (B6) 



p,q,r,s 



jT jM* 



By injecting the expression of kfj into [ifY (Eq lB4p . the intensity of the target peak is: 

[I^f = M*a^TS{xi-XT) (B8) 

where S denotes the Kronecker symbol. The muhiple scattering contribution (Eq IBSp results 
in an image of speckle whose mean intensity / [-^*^]"/ is given by: 



[in') = 2M^al (B9) 

The factor 2 comes from the fact that kf^ = kj^. The third and fourth terms, Iflj^* and 
jMjT*^ correspond to the interference between the signals associated with the target and 
the multiple scattering contribution. These signals are totally decorrelated, hence 

{If in = w'ln = (BIO) 

In average, the intensity of the echographic image, at the target arrival time, exhibits the 
following profile: a peak linked to the target at xi = Xt, of intensity M^cr|., buried in a 
speckle pattern whose mean intensity is 2M^o"|j-. 

If the maximum of the image, /max? is actually linked to the target, then its amplitude 
is given by: 

/max ^ E [/max] = J M^a'^ + 2M^al (Bl 1) 
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The quadratic mean of the image is: 



1 

1=1 



M3 (4 



20 



(B12) 



Upon normahzation (Eqj2^. /max is thus given by: 



'max 



M 



a. 



T 



(B13) 



At the hmit of detection, we can assume that M'^a^ » a\i » cr^. The vahdity of this 
approximation will be proved a posteriori by the final result. Eq JB13l simplifies into: 

~ Mar 



'max 



(B14) 



As to FB, the detection threshold corresponds to the condition /max > where a is given 
by Eql26] and depends on the PFA 7. Finally, we obtain the following detection criterion: 

aV2 



> 



(B15) 



CTM M 

This condition is reported in TabUl It indicates the signal-to-noise ratio ^ above which 
the main peak will correspond to the target (with a probability of false alarm 7) . Note that 
Eq JBlSI is only valid for a multiple scattering noise, which is spatially reciprocal. If we had 
dealt with an additional noise which does not respect this property, the detection criterion 
would be: 

0^1^ M 



2. The DORT method 



First, we deal with the classical DORT method, i.e when the SVD is applied directly to 
the raw matrix K°. If the first singular value A? is actually linked to the target echo, then 



its expected value is given by 



36|: 



E [A?] = Mar 



The quadratic mean of singular values is given by 36l |: 



M (4 + o 
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Upon normalization(Eq l22ll . the expected value of A? is thus given by: 
E 



The apphcation of the detection criterion A? > « leads to the following detection condition: 

— > ^ (B16) 



This condition is reported in TabUl If we compare it with the one obtained above for FB 
(Eq JBlSl ). we see that the DORT method is clearly more sensitive to noise than FB. The 
detection criterion varies as M~^^^ for the DORT method, whereas it is for FB. 

The argument is the same for the DORT method applied to K^, except that we have 
to take into account the action of the SSF. In Sec JIV Bl we have shown that the filtering 
of antidiagonals decreases the multiple scattering contribution by a factor \J^- Thus, the 
filtered matrix can be decomposed as: 

K^ = K^ + K^^ (B17) 

The matrix is linked with the residual contribution of multiple scattering. is a 

random Hankel matrix whose coefficients have a variance of 2a\j/M. The results obtained 
for K° can be applied directly to K^, taking into account the lower variance of multiple 
scattering signals. If the first singular value Af is associated to the target echo, then its 
expected value is still given by: 

E [Af ] = M(Jt 
The quadratic mean of singular values is given by : 



p=i 



2a 



Finally, upon normalization (Eq l22l) . the expected value of Af is given by: 



E<{ Af = . M ^ ^ , for M4 « alj 



cr|, ^ M ar 

The validity of the approximation Ma^ << a\j will be proved a posteriori by the final 
result. The application of the detection criterion Af > a leads to the following detection 
condition: 
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This condition is reported in TabUl If we compare it with the one obtained for FB (Eq JB15P 
and for the classical DORT method (Eq lB16() , we see that the SSF improves the detection 
condition by a factor \fM compared to the DORT method and hence, reaches the level 
of performance of FB. Actually, it is even slightly better than FB since the threshold a 
is inferior for DORT applied to K^. This approach succeeds in detecting the target for 
higher noise-to-signal ratios, compared to FB. This better performance is reinforced by the 
robustness of the SSF to aberration (see Sec JVlip . which has been neglected in this Appendix. 



APPENDIX C: EFFECT OF THE SSF ON ABERRATION 



This appendix deals with the effect of the SSF on aberration. More particularly, we 
want to express the coefficients kf.^ of the filtered matrix K^, when aberration exists. The 
coefficients of matrix K are given by (Eq l29l) : 



jk 



2R 



exp 



2Rr 



The coefficients dij form the distortion matrix D. The first step of the SSF consists in the 
rotation of data described in Sec lIV Al It results in two antidiagonal matrices Ai et A2 
(Eqsj6]&[71). Let A^ et A^ be the two antidiagonal matrices built from matrix D according 
the same process of Sec. IIV A[ 

Af = [af„,,] of dimension (2M - 1) x (2M - 1), 

such that af [n, v] = d[u + v - l,v - u ^ 2M - 1] (CI) 

A° = [ag^J of dimension (2M - 2) x (2M - 2), 

such that a2[u,v]=d[u + v,v-u + 2M -1] (C2) 

From now on, we will call indifferently A, the matrices Ai and A2, and A^, the matrices 
Af and A^. 

The coefficients of A can be expressed with the coefficients of A^: 









exp 



jk 



2R 



(C3) 



with 



Vu 



et y^, - 



72 



72 
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The next step of the filter consists in projecting the columns of A along the characteristic 
space of single scattering, generated by the vector S whose coordinates are: 



exp 



jk 



yl 



2R 



^-1/2 



The coordinates of the vector P, result of this projection (Eq fT2|) . can be expressed as : 



Pv 



u=l 

1 ^ 

^ y^exp 

' 1 ^ 













«™ exp 


2Rt 


exp 



jk 



2Rt 



,D 



exp 



jk 



2R 



The filtered matrix is finally obtained by multiplying the column vector S with the line 
vector P (EqJT6!). As a result, the coefficients a^^ are given by: 



F 



a 



1 ^ 



^ u=l 



exp 



jk 



yl 

2Rj 



exp 



jk 



{y, - V2Xt) 



21 



2R 



■T 



(C4) 



If we compare the expressions of auv (Eq ICSp and a^^ (Eq JC4p . we see that the SSF averages 
the distortion coefficients along each column of A. This average corresponds to the term 

Z Eti «™] in EqEa 
Afterwards, two matrices Af and A2 (containing the filtered antidiagonals) are obtained. 
The last step consists in constructing the filtered matrix K^, as described in Sec JIVCl : 



if (/ — m)/2 is an integer 

then, k^[l, m] = af [(/ - m)/2 + M, (/ + m) /2] 
if {I — m)/2 is not an integer 

then, k^[l, m] = [(/ - m - l)/2 + M, (/ + m - 1) /2] 



(C5) 
(C6) 



By injecting Eq JC4l into the two last equations and reversing the change of coordinates, we 
obtain: 



k [l,m] = ez+m-iexp 



jk- 



2R 



exp 



jk- 



{Xrri — Xj 

2Rt 



(C7) 
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where the coefficients e^+m-i are defined as: 

if (/ — m)/2 is an integer, 



^ 2M~1 

then, ei+m.^i = _ ^ ^ af [n, (/ + m)/2] 

u=l 

if (/ — m)/2 is not an integer, 

^ 2A/-2 

then, ei+™_i = _ ^ ^ [m, (/ + m - l)/2] 



(C8) 



(C9) 



We can finally express coefficients Cy as a function of distortion coefficients dij, using Eqs JCll 



if w is an odd number, 

2M-1 



2M - 1 ^ 



then, e„ 



u=l 

if y is an even number, 

^ 2A/-2 

then, e„ 



2M-2 



f — 1 V — 1 

u + ——,— M + 2M-1 



^ + |, ^ - « + 2M - 1 



(CIO) 



(Cll) 



Hence, the coefficients correspond to the average of distortion coefficients along each 
antidiagonal of the matrix D. The latter equations can be rewritten as: 



if V is an odd number 
then, Cy = ( d 



V — 1 V — 1 
u+ , U + 2M - 1 



?i=l,...,2M-l 



if V is an even number. 



then, e„ 



d 



V V 

n + -, u + 2M - 1 

, 2' 2 



u=l,...,2M-2 

where the symbol < . > denotes an average over the variable in the subsript. 



(C12) 



(C13) 



APPENDIX D: DETECTION CONDITION IN PRESENCE OF ABERRATION 



In Sec JVIIl we have already shown that the filtered can de decomposed as follows: 

MatrijToFrank 1 Random Hankcl matrix 
(K^) is the mean of K^, its coefficients are given by: 



(kLif)) = (e«+m-i(/))exp 



jk 



2Rt 



exp 



jk- 



2Rn 
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The norm of its entries is uniform and equal to the mean of the distortion coefficients, (e^). 
The matrix (K^) is of rank 1. 

AK^ corresponds to a perturbation(not necessary small) linked with the fluctuations of 
distortion coefficients Cp. Its coefficients are given by: 



^kimif) = [ez+m-i(/) - (ei+„_i(/))]exp 



jk- 



2Ri 



exp 



2R 



■T 



The coefficients Skj^{f) are random variables whose standard deviation is std [e^]. As seen 
previously, the matrix AK^ displays the same statistical behavior as a random Hankel 
matrix. 

The aim of this appendix is to determine the ratio | (cp) | / std [cp] , above which the S VD 
will succeed in extracting the matrix (K^) along the first singular space. To that aim, we 
will use RMT once again. 

We will rely on the same kind of argument as in Appendix IB 2[ Indeed, an analogy 
can be made with the detection condition found for the signal-to-noise ratio ax/uM- The 
signal amplitude ar corresponds here to the mean value of distortion coefficients |(ep)|. The 
standard deviation of Cp plays the role of "noise" : std [cp] -v^ cta/ . Note that the matrix 
we consider here is of size x "y- If the first singular space of corresponds actually to 
the "non distorted" matrix (K^), then the expected value of Af is given by (see Appendix 
El: 

E[Af] =|^Kep)| (Dl) 
The quadratic mean of singular values is given by (see Appendix IB 2^ : 

^ E [K] ' - (l(e,)r + var[ep]) (D2) 

Upon normalization (Eql22l). the expected value of Af is thus given by: 



The validity of this approximation |(ep)|^ << var [ep] will be proved a posteriori by the final 
result. As is a Hankel matrix, the detection criterion is Af > a. It leads to the final 
detection condition : 



>^= (D4) 



std [cp] 
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